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We report on recent progress in employing the Highly Improved Staggered Quark (HISQ) action 
introduced by the HPQCD/UKQCD collaboration in simulations with dynamical fermions. The 
HISQ action is an order a 2 Symanzik-improved action with further suppressed taste symmetry 
violations. The improvement in taste symmetry is achieved by introducing Fat7 smearing of the 
original gauge links and reunitarization (projection to an element of U(3) or SU(3)) followed 
by Asq-type smearing. Major challenges for calculating the fermion force are related to the 
reunitarization step. We present a preliminary study of the HISQ action on two 2+1+1 flavor 
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ensembles. 
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1. Introduction 

The Highly Improved Staggered Quark (HISQ) action developed in Ref. [|l|] is an 0(a 2 ) 
Symanzik-improved action for which additional suppression of taste-exchange interactions is 
achieved by replacing the original gauge links U in the Dirac operator by 



where intermediate sets of links V, W, X are defined as 

• J^i - smearing level 1 (Fat 7): V = ^\U, 

• - reunitarization: W = WV, 

• J^2 - smearing level 2 (Asq): X = j^W '. 

A new feature of the HISQ action compared with asqtad is the reunitarization step that gives 
an extra contribution to the fermion force. 

2. Fermion force and reunitarization 

In molecular dynamics simulations we are sampling an ensemble of gauge configurations 
weighted by exp(— S), where the action S = S g +Sf is split into gauge, S g , and fermionic, Sf, parts. 
As usual, the integration over Grassmann variables is performed and then a set of pseudo-fermion 
fields <t> is introduced, resulting in the fermionic part of the form 5/ ~ (<t>\(M 1 < (U)M (U))~ Nf / 4 \<t>) . 
One can consult the details of the algorithm in [Q] and recent ideas on efficient evaluation of the 
fermion force for HISQ in 

The fermion force is calculated by taking the derivative of the action 5/ with respect to funda- 
mental gauge links U using the chain rule along the lines of Refs. [^J, Schematically: 



The following parts are the same as for the asqtad action: dSf/dX, dX/dW, dV /dll, while the 
contribution from the reunitarization step dW /dV is new. 

We have experimented with projecting to SU(3) and U(3) groups and found that pion spectrum 
measurements with valence HISQ and sea asqtad quarks show no difference, as one would expect 
on physical grounds. Therefore in our dynamical HISQ simulations we always perform projection 
to the group U(3). While keeping the physics the same, this has two advantages: 

1. Different methods of projection, namely polar decomposition and trace maximization, give 
identical results. (The same is not true in the SU(3) case.) 

2. When projecting from U(3) to SU(3) by making the determinant of the matrix equal to 1, one 
needs to choose among three possibilities for the phase. In dynamical simulations one needs 
to uack this phase to make sure it changes smoothly for each link. Otherwise discontinuous 
changes lead to rapid changes in the action similar to those we describe below. 



(1.1) 



dS f _ dS f dX dW dV 
~dU ~ ~dX dW ~dV dU ' 



(2.1) 
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We adopt the following method of U(3) projection. For a complex matrix V the matrix H = 
V^V is Hermitian and W = VH^^ 2 is unitary. We first calculate H~ l l 2 by using the Cayley- 
Hamilton theorem in a manner similar to the approach in Refs. 

H-V 2 = f I+fiH+f 2 H 2 , (2.2) 

where are functions of c,- = Tr(H l+i )/(i + 1), i = 0, 1,2. To evaluate the derivative of H~ l l 2 
with respect to H one needs to know the derivative of fj with respect to H. This can be performed 
analytically by applying the chain rule and making use of the coefficients explicitly calculated in 
Ref. After dH~ l l 2 / dH is known we have trivially: 

dW dV dH~^ 2 dH 

= — H i +V (2 3) 

dV dV + dH dV { ' 

(We write Eq. ( j2~3| ) schematically, but in fact each matrix-matrix derivative is a rank 4 tensor with 
indices contracted in such a way that the resulting expression on the right hand side is a rank 
4 tensor again. An excellent review on how to deal with such derivatives can be found in [Q].) 
Thus the entire derivation can be performed analytically. (We also implemented a finite difference 
scheme and rational function approximation but found them much less accurate than the procedure 
described here.) 

To understand the possible behaviour of dW /dV qualitatively let us consider the U(l) case for 
the moment. Then V = re' 9 and W = V (V^V)^ 1 ^ 2 = e l9 , and the derivative 

dW_/dW\ _d(W,V Jf ) _d(W,V J{ ) d(r,0) 1 
~dV~\dv) vf ~d(V,V^)~ d(r,d) d(V,V^~2? ( ' } 

is large when r is small, as might happen when the links U are locally disordered, and the first 
stage of smearing results in a small smeared link V. For the matrix case the derivative is dominated 
by the smallest eigenvalue of V. To derive Eq. ( |2.4[ ) we applied the method of Jacobians and the 
notation is similar to the one used in Thermodynamics. 



3. Dynamical HISQ simulations 

We run the RHMC algorithm with five pseudofermion fields. The first pseudofermion im- 
plements the ratio of the determinants for two light and one strange quark to the determinant for 
three unphysical heavy quarks ("UHQ") with mass amuuQ = 0.2. The next three each implement 
the determinant for one UHQ, and the final pseudofermion implements a physical charm quark, 
including the mass correction to the Naik term to first order in m. 

At the time we did our studies, one loop fermion corrections to the gauge action were not yet 
known, so we used the coefficients in our gauge action appropriate for the asqtad fermion action. 
(The coefficients of the gauge action with 1-loop corrections due to HISQ have been calculated and 
are now available [§J.) 

Our integration algorithm is the "3G1F" ("three gauge steps, one fermion step") algorithm, 
with the Omelyan integrator used for both gauge and fermion forces. Our convention for the step 
size is that each application of the fermion force is one step. Note that since the Omelyan integrator 
shifts the time at which the force is calculated alternately forward and backward, the full period in 
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xN t 


Conf. 


P 


ami 


am s 


am c 


At 


a, fm 


20 3 x 64 


40 


6.75 


0.010 


0.050 


0.600 


0.04167 


0.127 


28 3 x 96 


50 


7.07 


0.007 


0.035 


0.420 


0.03125 


0.093 


48 3 x 144 


6 


7.47 


0.004 


0.020 


0.240 


0.01250 


0.060 



Table 1: Dynamical HISQ 2+1+1 ensembles. 




Figure 1: The histogram of the change in the action. 



simulation time of the integrator is twice what we call the step size. For example, when we run 
a fifty-step trajectory this means 25 cycles of the Omelyan integrator, involving 50 fermion force 
calculations and 150 gauge force calculations. The parameters of the three 2+1+1 flavor ensembles 
we have run are compiled in Table [I]. The lattice spacing is calculated by measuring the ratio r\ja 
from the static quark potential and using the value r x =0.318 fm 

While most trajectories ran smoothly, we found many trajectories with large jumps in the 
action as the integration proceeded. In Fig. [T] we show the histogram of the change in the action, 
plotted on a logarithmic scale, over several time units for the a = 0.127 fm ensemble. The long tail 
"outliers" indicate instantaneous jumps in the action, which we investigate further. 

Let us denote the norm of a matrix A by: 

When we calculate the fermion force at each time step as an anti-Hermitian matrix defined on each 
link, we evaluate its norm and find the maximum value over the lattice: | |F | \ max . Also, at each time 
step we calculate the determinant of V (Fat7 smeared) links and find the minimum value over the 
lattice: |detV|„„>,. Time histories of these two quantities are shown in Fig. ||. Large values of the 
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0.012 




Figure 2: The time history of | detV | m! „ (top) and | \F\ \ max (bottom) for the a = 0.127 fm ensemble during 5 
time units (= 5 x 24 = 120 time steps). 



fermion force accompany small values of the determinant (or small eigenvalues) of Fat7 smeared 
links. Thus, when calculating smeared links by adding different paths one may by chance produce 
a matrix V with an eigenvalue close to 0, which in turn leads to a large derivative dW / dV (as in the 
U(l) example considered earlier) that results in a large fermion force. Our integration algorithm 
has a finite step size, so it is not able to integrate such "spikes" in the force smoothly, leading to 
instantaneous jumps in the action that we see as "outliers" in the action histogram Fig. [j]. 

We investigated how this situation changes when we go to finer (smaller lattice spacing) en- 
sembles. Since configurations become smoother, spikes in the force become less severe, as seen in 
Fig-! 



4. Pion splittings on dynamical HISQ configurations 

The effect of suppression of the taste-exchange interactions in HISQ was investigated in [jl|] 
by measuring the pion spectrum with valence HISQ on sea asqtad configurations. Here we report 
on similar measurements performed on dynamical HISQ configurations for the first two ensembles 
of Table |. 

It is convenient to define a dimensionless quantity which is almost independent of quark mass: 

A = (Ml-M 2 G )ri, (4.1) 

where Mq corresponds to the Goldstone pion and M % refers to one of the other seven pion tastes in 
Tables |2| and ||. We calculate A for comparable asqtad and HISQ configurations, and then the ratio 

A ASQ 



R 



(4.2) 
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Figure 3: Fermion force ||F|| maA for different ensembles versus time in time units (1 time unit corresponds 
to 24 steps for a = 0. 127, 32 steps for a = 0.093 and 80 steps for a = 0.060 ensemble). 
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I 


0.359(5) 
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0.079(4) 


0.0280(61) 


2.8(6) 



Table 2: Pion spectrum on a = 0.127 fm HISQ ensemble. 



shows how much the splittings decrease when going from asqtad to HISQ. The values of R for 
different pion tastes are shown in the last column of Tables |2| and [3|. The statistical errors are rather 
large since for HISQ we have about an order of magnitude fewer configurations than for asqtad. 
The number of configurations used for measurements is indicated in parentheses in the headers of 
the second and third column. The overall trend is however clear and in agreement with Ref. [jl|] : 
about a factor of three improvement in taste symmetry for the HISQ action relative to asqtad. 



5. Conclusions 



In dynamical HISQ simulations with typical parameters, we found that smearing may produce 
a smeared link V with a small eigenvalue that dominates the derivative of the reunitarized link W, 
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Pion taste 
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0.2256(11) 
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0.0028(4) 


3.0(4) 


Yi 


0.2311(15) 
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0.0031(5) 


3.5(6) 


7o 


0.2318(10) 


0.1509(11) 


0.0109(5) 


0.0038(4) 


2.9(3) 


I 


0.2398(25) 


0.1522(27) 


0.015(1) 


0.0042(9) 


3.5(8) 



Table 3: Pion spectrum on a = 0.093 fm HISQ ensemble. 

dW /dV, and gives a large contribution to the fermion force. Such "spikes" in the force integrated 
with finite time steps give "jumps" in the action that decrease the acceptance rate of the RHMC 
algorithm. This problem was noted in Ref. ^ and is probably related to topological defects, "dis- 
locations" that manifest themselves in plaquettes with low values lying on the tails of the plaquette 
distribution. We found that going to finer ensembles (smoother gauge configurations) reduces the 
number of spikes and partially cures the problem. 

On configurations generated with the HISQ action for dynamical quarks we measured the 
staggered pion spectrum and found that pion splittings decrease by a factor of three or more, con- 
firming the result of Ref. [|l|], where pioneering tests were done with valence HISQ fermions on 
configurations with asqtad sea quarks. 
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